Analysis of parametric biological models with non-linear 

dynamics 



In this paper we present recent results on parametric analysis of biological models. The underlying 
method is based on the algorithms for computing trajectory sets of hybrid systems with polynomial 
dynamics. The method is then applied to two case studies of biological systems: one is a cardiac 
cell model for studying the conditions for cardiac abnormalities, and the second is a model of insect 
nest-site choice. 

1 Introduction 

Computer aided modelling and analysis have been intensively used to study biological systems. Different 
models have been combined in a unified framework to capture the complex behavior of biological phe- 
nomena. Recently, the hybrid systems paradigm has been applied to various biochemical networks with 
switching behavior. On the other hand, the main difficult in biological model validation is the precision 
of predictions provided by the model, in comparison with experimental data. Quantitative information 
on kinetic parameters is thus hard to obtain and, in addition, experimental data are often of qualitative 
nature, which poses a major problem to traditional numerical methods. In this work, we try to address 
this problem by extending the numerical approach to systems with uncertain parameters. 

More concretely, we exploit the ability of handling uncertainty of the reachability computation tech- 
niques for hybrid systems, a mathematical model for describing the interaction between a discrete and a 
continuous process. Roughly speaking, the goal of reachability analysis is to study the set of all possible 
trajectories of a dynamical system. Well-established results for piecewise linear systems are available 
and they have been successfully applied to numerous biological systems (for example (2l [12j |9l [T5l0 . 
Nevertheless, nonlinear systems still remain a challenge. 

In particular we focus the following computation problem: given a set of initial states in W, compute 
the trajectory set of a parametric discrete-time dynamical system described by the following difference 
equation: x(k + 1) = n(x(k),p) where % : W x W" — > W is a multivariate polynomial. This problem can 
be seen as an extension of numerical integration, that is, solving the equations of the dynamics with sets, 
that is x(k) and x(k + 1) in the equation are subsets of W (while they are points if we only need a single 
solution, as in traditional numerical integration). 

Our interest in polynomial systems is motivated by their applicability in modeling a variety of phenom- 
ena in bio-chemical networks. To handle continuous-time models, a time discretization can be needed 
and conservativeness of approximation needs to be guaranteed (which is the topic of our current re- 
search). It is however important to note that discrete-time systems can also be directly used to model 
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many biological systems, since experimental data are often obtained by sampling the biochemical reac- 
tion outputs, and in addition, such models can be readily used for computer aided analysis and numerical 
simulation. 

This problem for polynomial systems was previously considered in the work ll24l [Toll , which was 
inspired by the techniques from Computer Aided Geometric Design (CADG). In this paper, we pursue 
the direction by exploiting further the special properties of a technique from CADG, namely the Bezier 
simplex and Bernstein expansion, we only need to solve linear programming (LP) problems instead 
of polynomial optimization problems. In this paper, with a view to biological models, we propose an 
extension of these techniques to systems with uncertain parameters. We also describe a new techniques 
specialized for multi-affine systems 

The paper is organized as follows. In Section [2] we introduce basic definitions of set integration. We 
then focus on the problem of computing the image of a set by a parametric polynomial and propose 
two methods: one is an extension to parametric systems of the method using the Bernstein expansion, 
and the other method is based on a property of multi-affine functions. Some experimental results on two 
biological systems are reported in Section[5} one is a model of cardiac abnormalities, the other is a model 
of insect nest-site choice. 

2 Set integration 

Let R denote the set of reals. Throughout the paper, vectors are often written using bold letters. Excep- 
tionally, scalar elements of multi-indices, introduced later, are written using bold letters. Given a vector 
x, %i denotes its i th component. 

We use SS to denote the unit box anchored at the origin, that is = [0, 1]". We use % to denote a vector 
of n functions such that for all i € {1, ... ,n}, 7T, is an rc-variate polynomial of the form 7T; : W — > M.. 

We consider a parametric discrete-time dynamical system: 

X(*+1) = JE(X(*),P) (1) 

where x £ W is the state variables, pGPC W" is the vector of uncertain parameters. The set P is called 
the parameter set. We assume that P is a convex polyhedron. The initial state x(0) is inside some set 
Xo C W, and Xo is called the initial set. 

Given a set X C W, the image of X by % with all possible values of the parameters, denoted by 7r(X), 
is defined as follows: n(X) = {(7Fi(x,p),. . . ,7T„(x,p)) | x G X,p G P}. 

When starting from a set of initial states, the dynamical system ([T]) generates a set of solutions, and at 
each step the set of all visited states is the image by 71 of the previous set. We thus focus on the problem 
of computing the image of a polyhedron X by the polynomial %. 

One way to charaterize this set of solutions is to use special convex polyhedra with fixed geometric 
form, called template polyhedra ll23l ITI. Indeed, by choosing the templates one can control both the 
precision and the geometric complexity of the approximations. A template is a set of linear functions 
defined by an / x n matrix H and a real- valued vector c E Mf. A template polyhedron is then defined by 
considering the conjunction of the linear inequalities of the form (H,c) = {x | f\ j=l iH'x < c,}. The 
vector c is called a polyhedral coefficient vector. 

The template matrix H is assumed to be given; to over-approximate the trajectory set, the polyhedral 
coefficient vector c6K' must satisfy 

n(X)Q{H,c). (2) 



18 



Analysis of parametric biological models with non-linear dynamics 



To determine the coefficients c = (c\ , . . . , c n ) T so that the template polyhedron (H, c) over-approximate 
%(X), we can formulate the following optimization problems: 

Vi € {1,. . . ,l}, Ci = max(£2 =1 #!7t*( X )) (3) 
subj. to x £ I. (4) 

where H' is the i row of the matrix H and H' k is its k th element. Note that the above functions to 
optimize are polynomials. This polynomial optimization problem is computationally difficult (see for 
example [13] and references therein). We propose two alternative solutions: 

• Use a linear relaxation of the above optimization problems, so that we can take advantage of well- 
developed linear programming techniques ll22ll . Indeed, the Bernstein expansion can be used to 
compute affme bound functions of polynomials, as shown in the next section. 

• If 71 is multi-affine, we can exploit its particular properties that allows using corner analysis to 
reconstruct the set of all trajectories. 

In the following section, we describe two methods for approximating the image of a box domain by a 
polynomial, and we then show how to extend these methods to polyhedral domains. 



3 Approximating the image of a box domain 

3.1 Using the Bernstein expansion 

The essence of this method, proposed in lfl4l . for computing an affine bound function is to find a hy- 
perplane that is close to all the control points, using linear least squares approximation. In this work, 
we extend this method to parametric polynomial functions. Before presenting the Bernstein expansion, 
we remark that computing convex lower bound functions for polynomials is a problem of great interest, 
especially in global optimization. The reader is also referred to ll20ll for more detail on the approach 
based on the Bernstein expansion. 

To discuss the Bernstein expansion of polynomials, we use multi-indices of the form i = (h,i2, ■ ■ ■ ,i n ) 
where each i ; - is a non-negative integer. Given two multi-indices i and d, we write i < d if for all 

j € {1, . .. ,n}, ij < dj. Also, we write j for (ii/di,i2/d2, .. . ,i n /d„) and \ for the product of binomial 

coefficients (V I ( !M • • • ( H I ■ 
\dj \d 2 J \d n J 

The parametric polynomial 71 : W x W — > W can be represented using the power base as follows: 

^( X >P) = 52 a i(P) x> 
ie/ d 

where a, is a function R m — > R; i and d are two multi-indices of size n such that i < d; la is the set of all 
multi-indices i < d, that is 1^ = {i | i < d}; x' represents the monomial x{x^ ■ • The multi-index d is 
called the degree of %. 

The polynomial % can also be represented using the Bernstein expansion. For x = (jci , . . . ,x„) S W, 
the \ tb Bernstein polynomial of degree d is defined as follows: 



^d,i(x)=j3 dl , il (x 1 )... j 8 d „, i „(x„) 
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where for a real number y, fi^^iy) = —y Aj '■')• 

Then, for all x G S3 = [0, 1]", the polynomial % can be written using the Bernstein expansion as follows: 

7 r(x) = £b i (p)^ d , i (x) 

ie/ d 

where for each ie/d the Bernstein coefficient bi(p) is defined as: 

bi(p) = Llaj(p). (5) 
j<i IjJ 

With respect to our problem of computing the image of a set by a polynomial, the Bernstein expansion 
is of particular interest since their coefficients allow to derive useful geometric properties. We now show 
that these coefficients can be used to efficiently compute an affine approximation of the polynomial %. 

We denote by = {V | 1 < j < nt,} be the set of all the multi-indices, tit, is thus their number. The 
set of all control points are denoted similarly. Let A be a matrix of size n/, x [n + 1 ) (n is the number of 
state variables of the dynamical systems in question) such that its elements are defined as follows. For 
all 1 < j < nb an d 1 < k < n, 



andA „+i = 1 - 

We pick the centroid point p c in the parameter set P. Let £ be the solution of the following linear least 
squares approximation problem: 

A T AZ=A T b(p c ). 

Then, the affine function 

n 
k=l 

corresponds to the "median" axis of the convex hull of all the control points. It thus suffices to shift it 
downward by the amount: 

- i j 

5 = ffw{l( 7 )-b J p I < j < n b A p 6 P\. (6) 
d 

If the function % is linear in the parameters p, then the above optimization problem is a set of linear 
programs and can thus be solved efficiently. This results in a lower bound function 

l(x) = I(x)-8, for all xe^. 

We now show how the above affine bound functions can be used to solve the optimization problems ([3]) 
in order to determine the coefficients of a template polyhedron over- approximating the reachable set. The 
functions to optimize in Q can be seen as the compositions of polynomials Tfy. Since every coefficient 
H' k is constant, each 

a ; '(x)=£Li^« 

is simply a polynomial and we can compute its bound functions. The template polyhedral coefficients 
can hence be computed by solving the following optimization problems: 



Vi G {!,...,/}, Ci = max(o'(x)) subj. to x G X. 



(7) 



20 



Analysis of parametric biological models with non-linear dynamics 



Linear least squares approximation 



Translation 





Figure 1: Computation of the affine lower bound function using linear least square approximation; the 
vertical line segments represent the intervals of the Bernstein coefficient values. 



Example. We illustrate the method for computing the bound functions with an example. We consider 
a one dimensional polynomial of degree 5: 

7l(x) =p(l-x + 2x 4 ) + 3x 2 -x 3 -2.5x 5 

where p is a scalar parameter included in [0.5, 1.5]. Using the equation Bernstein coefficients are the 
following: 

b(p) = (p,0.8p,0.6p + 0.3,0.4p + 0.8,0.6p+1.4,2p-0.5). 
Next we compute the matrix A and the vector b(p c ) where the centroid p c is equal to 1: 

0.2 0.4 0.6 0.8 1 



11 1 1 1 1 J ' 

b(p c )=b(l) = (1,0.8,0.9,1.2,2,1.5), 
which gives us the following linear system: 



2.2 3 \ / Ci 
3 6 U 



4.34 
7.40 1 ' 



By solving this linear equality system, we obtain following least square approximation 

f(x) = C1X + C2 = 0.9 143x + 0.7762. 
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Then it follows from Q that 8 = 1.1905 and we get the affine lower bound 

l(x) = f(x) -8 = 0.9143x - 0.4143. 

The left plot in Figure [T] shows the linear least square approximation computed for the polynomial 
in question. The other one represents the final affine function computed after shifting downward the 
previous affine function. 

3.2 Using the property of multi-affine functions 

In this section, we propose a method specialized for multi-affine systems, a particular case of polynomial 
systems. Multi-affine systems are systems composed by polynomials which are affine by each of theirs 
variables, i.e. if d is the degree of an affine system 71, dk < 1 for all k G {1, . . . We will exploit the 
following property of such systems to achieve better computational efficiency. 

Theorem 1 Given a hyper-rectangle X C W, let Vx be the set of its vertices. If K is a multi-affine 
function, then 

k{X) C conv{n{y) | v G V x } 

To prove the theorem, we first prove an intermediate result. 
We denote the line segment connecting two points x and x' by s(x,x'). 

Lemma 1 If two points x and x' that differ only in one coordinate ( that is there exists only one axis 
i € { 1 , . . . , n} such that Xj ^ x' i and for all other axes j ^ i Xj = x'j), then for any point y in the line 
segment s(x,x'), its image ?r(y) can be written as a linear combination ofn(x') and 7l(x). 

Proof. It is easy to see that any two points in s(x,x') also differ only in the axis i, (that is s(x,x') is 
parallel to the axis xi). The i coordinate of any point y 6 s(x,x') can be expressed as y,- = Aix, + A2X- 
where X\ , A2 > and X\ + = 1. In addition, if the domain of the function % is restricted to the line 
segment s(x,x'), only the variable xi varies while all the variables remain unchanged and the function % 
hence becomes linear in x\. It then follows that 7r(y) can be written as a linear combination of 7t(x') and 

7T(x). ■ 

We proceed with Theorem [T] Let x is a point inside X. Let y and z be the projections of x on two 
opposite faces of X perpendicular to the axis x n . It follows from Lemma[T]that n(x) can be written as 
a linear combination of 7i(y) and 7i(z). Furthermore, y can now be treated as a point inside a hyper- 
rectangle in (n — 1) dimensions (jci , . . . ,x ni ), and by successive projections until dimension x\ we can 
prove that 7i(x) can be written as a linear combination of the vertices of X. ■ 
Another proof of this well-known property of multi-affine functions can be found in [4J 
Note that the above theorem is only applicable to the sets which are axis-aligned hyper-rectangles. 
Hence, even if the initial set satisfies this condition, after the application of %, the resulting set is a 
general polytope and we need to approximate it by an axis-aligned hyper-rectangle. Using the theorem, 
to compute 7t(X) it suffices to compute the images of the vertices of X and then take the convex hull of 
the resulting points. 

Before continuing, we remark that a number of different methods for the multi-affine systems have 
also been developed in J3j|5l. These methods are however based on a rectangular partition of the state 
space, while we allow reachable sets to be represented by unions of polytopes. 
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4 Approximating the image of a polyhedral domain 

As mentioned earlier, the above described methods can be applied only when the set X is inside the unit 
box 38 anchored at the origin. To extend it to polyhedral domains, we transform the polyhedra to the 
unit box by two methods: (1) via an (oriented) box approximation, and (2) by rewriting the polynomials 
using a change of variables. 

If we over- approximate X with a box B, it is then possible to derive a formula expressing the Bernstein 
coefficients of 7t over B. However, this formula is complex and its representation and evaluation can 
become expensive. 

We alternatively consider the composition of the polynomial % with an affine transformation z that 
maps the unit box to B. The functions resulting from this composition are still polynomials, for which 
we can compute their bound functions over the unit box, using the formula Q of the Bernstein expansion. 
This is explained more formally in the following. 

Let B be the bounding box of the polyhedron X, that is, the smallest box that includes X. The affine 
function X that maps the unit box 38 to B can be easily defined as: t(x) = dictg(X)x + g where gGl" 
such that gi = lj, and diag(X) isanx/i diagonal matrix with the elements on the diagonal defined as 
follows: for each i G {1, . . . ,n}, A,- = hi — lj. The composition y = (n o t) is defined as y(x) = 7t(t(x)). 
Note that 7= nor. Then, n(X) C y(38). 

Since multi-affine functions are not closed under the above composition. They can become quadratic. 
In this case, we use the following blossoming principle to transform a polynomial to an affine function. 

Theorem 2 (Blossoming principle) For any polynomial % : W — > W of degree d, there is a unique 
symmetric multi-affine map w : W x . . . x M." — > M m such that for all x£l" 

w(x,...,x) = 7r(x). 

The proof of the above theorem is standard and can be found in |[T6ll . We can thus compute the image 
7t(B) by the convex hull of all the images of the vertices of X by the multi-affine function w. Nevertheless, 
it is important to note that when x inside an rectangle, the arguments of the blossom varies only on the 
diagonal of the corresponding rectangle in the augmented space M. nd . Thus the convex hull is an over- 
approximation. 

5 Experimental results 

We implemented the above described algorithms, and their time efficiency was also evaluated by consid- 
ering a number of randomly generated polynomials, which shows that the algorithms can handle general 
polynomial systems in up to 11 dimensions with reasonable computation time. In the following, we 
demonstrate their application to two biological systems. 

5.1 A model of insect nest-site choice 

We study a model of a decision making mechanism used by a swarm of honeybees to select one among 
many different nest-sites [6]. This is built upon classical mathematical models of the spread of diseases, 
information and beliefs. The bee population is divided into 5 groups: 

• X: neutral bees that have not chosen a site 
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• Y[ : evangelic bees dancing for the site 1 

• F 2 : evangelic bees dancing for the site 2 

• Z\ : non-evangelical bees which have been converted to the site 1 

• Z 2 : non-evangelical bees which have been converted to the site 2 

The bees dance not only to advertise the quality of a site but also to transmit to the other bees information 
about the direction and distance to the site. Non-evangelical bees stay idle but may take up dancing when 
they are stimulated by a dancing bee. For simplicity, we also use the names of the groups to denote the 
respective numbers of individuals in each group. 
The equations describing the evolutions of the variables are as follows: 

X(t) = -PiX(t)Yi(t)-P3K(t)Y 2 (t), 

= j 8iX(0Fi(0-7^i(0 + 5j8iFi(0Zi(0 + aj8iFi(0Z 2 (0, 
Y 2 (t) = faX(t)Y 2 (t)-YY 2 (t) + 8fcY 2 (t)Z2{t) + ap 2 Y 2 (t)Z l (t), 

2i(f) = 7Fi(0-5j8iF 1 (0Z 1 (0-ai82F 2 (0Z 1 (0, 
Z 2 (t) = yY 2 {t)-8p l Y 2 {t)Z 2 {t)-ap 2 Yi{t)Z 2 {t), 

where j3i and j3 2 are scalar parameters representing a measure of how vigorously the bees dance for the 
sites 1 and 2 respectively; a is the parameter representing the proportionality of switching back to the 
neutral state, and 7 is the proportionality of conversion from the dancing state to the idle state. The 
proportionality of conversion from the neutral state to any site F; is 1 and from the idle state to F, is 8. 

In this case study, we want to study the influence of the parameter /3 2 on the possibility of achieving a 
consensus on nest-site choice. 
Using the Euler discretization method we obtains the following discrete-time dynamics. 

X(k + 1) = X(fe)+ft(-/3 1 X(fc)F 1 (fc)-/3 2 X(fc)F 2 (fe)), 

Yi(k + 1) = Y\ (k) + h{frX(k)Y, (k) - 7F (*) + SftFj [k)Z x (k) + afaY Y (k)Z 2 (k)), 

Y 2 (k + l) = Y 2 (k) + h(p 2 X(k)Y 2 (k) - Y Y 2 (k) + 8$ 2 Y 2 {k)Z 2 {k) + afcF 2 (*)Z 1 (*)), 

Zi(* + 1) = Z 1 (^)+/ J (7F 1 (^)-5 j S 1 F 1 Z 1 -a j 8 2 F 2 Z 1 ), 

Z 2 (k + \) = Z 2 (k)+h(yY 2 (k)-8p l Y 2 (k)Z 2 (k)-ap 2 Y l (k)Z 2 (k)), 

We choose a discretization time h = 0.01. In the following, = 1.0, the initial set is a small box 
centered at (N,0, 0,0,0), where Af is the total number of honeybees arbitrary fixed to 1000. 

Figure [2] and Figure [3] show the evolution of the proportions of converted honeybees ((F,- + Zi)/N) 
during 6000 iterations, with different interval values of the parameter j3 2 . The blue sets (in lighter color) 
correspond to the proportions of bees supporting the site 1 and the red sets (in darker color) correspond 
to those supporting the site 2. These results were obtained using the method for multi-affine systems 
which took about 7.5s 

First we consider that at the beginning no honeybees dance for the site 2, in other words j8 2 = 0. After 
300 iterations the site 2 is discovered and its propaganda begins. We can observe two main comport- 
ments: a consensus and the absence of consensus. We consider that there is a consensus if the distance 
between the proportions of honeybees supporting the site 1 or 2 is large and tends to increase. 

Figure|2]shows the analysis using the parameters a = 0.7, 7 = 0.3 and j8 2 G [1, 1.2]. We observe that j8 2 
is marginally superior to j8i, however we can see that there is no consensus between the two evangelical 
groups. 
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Figure 2: No consensus is observed with a = 0.7, y = 0.3 and fa 6 [1, 1.2]. 

Figure [3] shows that a consensus can be observed if the measure of how vigorously the bees dance for 
sites 2 is increased. In this analysis the parameters stay the same except for fa £ [1.5,2.0]. We observe 
a clear consensus for the site 2 despite its lateness. 




time 



Figure 3: A consensus for the site 2 is observed with a = 0.7, 7= 0.3 and fa G [1.5,2]. 

Another way to obtain a consensus without modifying fa is to reduce a which corresponds to the 
proportionality of honeybees reconversion to another site. When a = 0.2, the late propaganda of the site 
2 leads to a consensus to choose the site 1 as shown in Figure |4] 

Finally, to compare the two proposed methods, Figure [5] shows the projection of the reachable sets 
on the variables Y\ and Y2 during the first 150 iterations, with a = 0.7, 7= 0.3 and fa = 1.2. We plot 
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time 

Figure 4: A consensus for the site 1 is observed with a = 0.2, 7= 0.3 and fa G [1, 1.2]. 

the results obtained by the two methods and we can see that the sets computed using the method for 
multi-affine systems are larger than those computed using the Bernstein expansion. We observe a gain 
of precision with the Bernstein method, however the computation time of this method is 25. 06s and is 
more than 250 times superior than the computation time of the other method which is 0.09s. 

5.2 A model of cardiac abnormalities 

The second case study comes from a simplified model of cardiac cells [ 15 ] which describes the condi- 
tions under which cardiac abnormalities (such as, ventricular tachycardia and fibrillation) can arise. The 
behavior of cardiac cells can be modelled using a Detailed Ionic Model (DIM) which are differential 
equation models of reaction-diffusion type. Such models are complex with a lot of parameters. An alter- 
native model is called Minimal Resistor Model (MRM), the parameters of which can be identified from 
either experimental data or from DIM-based simulation results. The main advantage of the MRM model 
is that it exhibits more accurate mesoescopic behavior and additionally allows faster simulations. From 
the MRM model, the authors of [ 15 ] derived an equivalent Minimal Conductance Model (MCM) which 
is exactly in form of a genetic regulatory network model (GRM). By approximating slow sigmoidal 
switches in the MCM model with a succession of ramps, a Hybrid Automaton ( ifTTl ) is obtained. 

In this case study, we consider the Hybrid Automaton model and the following problem: find the 
parameter ranges which correspond to the behaviors where the cardiac cells lose their excitability. These 
ranges can then be used to infer the corresponding parameter ranges in the DIM model. The hybrid 
automaton has two modes, denoted by loci and loci. 

• loci: 



du 

It =e ~ 81 
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Figure 5: Reachability analysis during 150 iterations using the two reachability methods, the black 
outlines set are computed with the Bernstein expansion method, the red outlines set correspond to the 
result using the method for multi-affine systems. 



The transition guard from loci to loci is u > 0.06. In this model, an action potential is a change in 
cell transmembrane potential u as a response to an external stimulus (current) e defined by e = 0.66 if 
t < 0.25 and e = if t > 0.25. The initial state corresponds to u = 0. In the above equations, g\ and 
g2 are parameters such that 1 < gi < 180 and 1 < gj < 10. The safety property we are interested in is 
stated as follows: the system, when at location loci, will never reach a state where u > 0.13. To this 
end, we add a fictitious location loci and let the transition guard from loci to loc3 be u > 0. 13. We also 
consider the two parameters as two additional variables with derivatives equal to 0. The problem is now 
to compute a set of possible values of g\ and g2 which guarantee that the system does not reach loci. 

The model hence becomes a Hybrid Automaton with 4 variables and their dynamics are described by 
multi-affine differential equations, which we discretize in time. Using backward reachability analysis 
starting from the unsafe set, we can over- approximate the set of (gi,g2) for which the system can reach 
the guard of the transition to loci. The figure [6] shows the projection on the variables u, g\, g2 of the 
reachable set, starting from the guard leading to the location loci. Its intersection with the hyperplane 
defined by u = (in grey) corresponds to the set of unsafe values of g\ and g2- This computation uses 
dynamical templates and static boxes, and the computation time is 1.525. The complement of this set is 
the set of safe parameter values. 



• loc2: 



du 

dt 
dv 

dt 



e-giu 



~Vg2\ 
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Figure 6: Backward reachability analysis, the change of color represents the change of location. 

6 Related work 

Our reachability analysis approach is similar to a number of existing ones for continuous and hybrid 
systems in the use of linear approximation. Its novelty resides in the efficient way of computing linear 
approximations. Indeed, a common method to approximate a non-linear function by a piecewise linear 
one, as in the hybridization approach dELLl for hybrid systems, requires non-linear optimization. 

Besides constrained global optimization, other important applications of the Bernstein expansion in- 
clude various control problems [19] (in particular, robust control). The approximation of the range of a 
multivariate polynomial over a box and a polyhedron is also used in program analysis and optimization 
(for example lPT8l |8l). In the hybrid systems verification, polynomial optimization is used to compute 
barrier certificates lT2"D . 



7 Conclusion 

The reachability computation methods we proposed in this paper combine the ideas from optimization 
and the properties of polynomials. These results can be readily applicable to hybrid systems with poly- 
nomial continuous dynamics. We are currently working on the problem of parameter synthesis using the 
analysis techniques from robust control. The next step will be a more intensive evaluation of the methods 
on biological systems. 
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Figure 7: Over-approximation of the parameter values (gi,g2) that lead to the unsafe location. 
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